%**************************************************************************
%
%           "Endogenous Liquidity and Capital Reallocation"
%      Cui, W., Wright, R., & Zhu, Y. (2024), Journal of Political Economy
%                       Last Modified: Feb 2024.  
%
%**************************************************************************

%%%% Run business-cycle simulations for the main (CWZ) model WITHOUT money

clc; clear; close all; format short


%%%% If 1, calibrate c_sig given c_mu, which is the version reported in
%%%% Table 8 of the paper; 
%%%% If 0 calibrate c_mu given c_sig
cali_c_sig = 1; 



%% 1. New calibration

disp('************> Model 6 (CWZ) Simulation without Money <***************')

load param_calibrated.mat
par                 = par_calibrated;

par.ETA             = 0.6168;   % recalibrate such that I/Y ratio is 20%; NOTE: this also means that sig should be recalibrated
par.log_eps_mu      = 0;        % log-normal capital productivity mean
par.log_eps_sig     = 0.5 / (1 - par.ETA); % log-normal capital productivity std.dev
par.S_PARAM         = par.log_eps_sig;
par.M_PARAM         = par.log_eps_mu;

fun  = define_function_fun(par);

disp('****** Calibration for the model case with credit only')
cali_step   = 1;
x0          = [0.08, 0.12, 1, 2];
x           = fsolve(@(x) eqm_ss_iid_cases_2_to_5(x, par, fun, 0, cali_step), x0);

[~, eqm] = eqm_ss_iid_cases_2_to_5(x, par, fun, 0, cali_step);

B                   = x(1);
par.ALPHA_0         = x(2);
par.CHI_q           = x(3);
par.XI              = x(4);
par.G               = eqm.G;

aqc_spending        = eqm.E_a_K  / par.ALPHA_0 / eqm.prob_A;
disp('Investment to output ratio')
disp(par.DELTA * eqm.K / eqm.Y * (1 - par.IM_share));

if cali_c_sig==1
    % calibrate c_sig given c_mu
    par.c_sig           = (log(eqm.surplus_net) - par.c_mu) / erfinv(2 * par.ALPHA_0 - 1) / sqrt(2);
else
    % calibrate c_mu given c_sig
    par.c_mu            = log(eqm.surplus_net)-par.c_sig  *sqrt(2)*erfinv(2 * par.ALPHA_0 - 1);
end



%% 2. Simulation

%%%% Related to dynamics
par.RHO_CREDIT      = 0.83;
par.RHO_A           = 0.83;
par.RHO_MU          = 0.83;
par.RHO_G           = 0.83;

x0 =  [0.08, 0.12];

x                   = fsolve(@(x) eqm_ss_iid_cases_2_to_5(x, par, fun, 0, 0), x0);
[~, eqm]            = eqm_ss_iid_cases_2_to_5(x, par, fun, 0, 0);


K_bar       = eqm.K_bar_K * eqm.K;
E_a         = eqm.E_a_K * eqm.K;
E_p         = eqm.E_p_K * eqm.K;
R_ratio     = (E_a + E_p) / (E_a + E_p + par.DELTA * eqm.K);
P_share     = E_p / (E_a + E_p);


%%%% store the steady state values
par.B_SS         = eqm.B;
par.Z_SS         = eqm.Z;
par.K_SS         = eqm.K;
par.H_SS         = eqm.H;
par.W_SS         = (eqm.B / (1 - par.ETA))^((par.ETA - 1) / par.ETA) * par.ETA;
par.L_SS         = max(par.EPS_L_BOUND, (eqm.Z/eqm.K - (1 - par.DELTA) * (1 - par.CHI_q - par.CHI_k)) / (1 - par.tau_k) / eqm.B);
par.Y_SS         = eqm.Y;
par.C_STAR_SS    = eqm.c_star;
par.K_BAR_SS     = K_bar;
par.q_capital_SS = eqm.q_capital;
par.q_money_SS   = eqm.q_money;
par.q_profit_SS  = eqm.q_profit;
par.surplus_SS   = eqm.surplus_net;
par.q_liq_SS     = eqm.q_liq;
par.E_a_SS       = E_a;
par.E_p_SS       = E_p;
par.R_ratio_SS   = R_ratio;
par.P_share_SS   = P_share;
par.DM_price_SS  = eqm.average_2nd_price;
par.prod_disp_SS = eqm.prod_std;
par.ALPHA_0_SS   = eqm.ALPHA_0;


cd dynamics_iid        
    disp('Simulatino with Credit and Productivity Shocks')
    if cali_c_sig==1
        par.SIGMA_CREDIT    = 0.817/100 * 1;
        par.SIGMA_A         = 3.297/100; % for productivity shocks only, multiply this by 1.21 to hit the volatility output
    else
        % use the following two lines if calibrate c_mu instead of c_sig
        par.SIGMA_CREDIT    = 0.797/100 * .584;
        par.SIGMA_A         = 3.286/100;% for productivity shocks only, set this to 3.286/100* 1.23
    end
    par.SIGMA_MU        = 0;
    par.SIGMA_G         = 0;
    save params par   
    dynare model_5_IRFs.mod nostrict
    oo_credit_prod = oo_;

    disp('Simulatino with only Productivity Shocks')
    par.SIGMA_CREDIT    = 0;
    if cali_c_sig==1
       par.SIGMA_A         = 4.02/100;
    else
       % use the following line if calibrate c_mu instead of c_sig
       par.SIGMA_A         = 3.286/100* 1.24;
    end
    par.SIGMA_MU        = 0;
    par.SIGMA_G         = 0;
    save params par   
    dynare model_5_IRFs.mod nostrict
    oo_prod = oo_;        
cd ..

disp('**********************> Results <**********************')


fprintf('output volatility with only productivity shock is %.4f \n \n', sqrt(oo_prod.var(1)))
fprintf('output volatility with both shocks is %.4f \n \n', sqrt(oo_credit_prod.var(1)))

std_dev_credit_prod = sqrt (diag(oo_credit_prod.var)) ./ sqrt(oo_credit_prod.var(1)); %the first position is output. Consumption starts from the 2nd
std_dev_prod = sqrt (diag(oo_prod.var)) ./ sqrt(oo_prod.var(1));

corr_w_output_credit_prod = oo_credit_prod.contemporaneous_correlation(1,:); 
corr_w_output_prod = oo_prod.contemporaneous_correlation(1,:); 


disp('-----------------------------')
disp('Table 8 result and more (last three rows, some are discussed in the paper):')
disp('Relative Std.Dev. and Correlation')
Table_8 = array2table([std_dev_prod(2:11), std_dev_credit_prod(2:11),...
   corr_w_output_prod(2:11)',  corr_w_output_credit_prod(2:11)'],...
   'RowNames', {'Consumption','Investment','Hours','TFP','R share', 'P share',...
    'Inflation','ALPHA_0','Prod disp','2nd Mkt Price'},...
   'VariableNames', {'std.dev. (A)', 'std.dev. (A&Credit)',...
    'corr (A)', 'corr (A&Credit)'})

